version 16.1
clear all
pause on

cap log close
log using cf_figures_xs.log, replace

global RESULTS ../output
global TEMP ../temp
cap mkdir ../temp/results_exhibits_xs
cap mkdir ../output/results_exhibits_xs
cap mkdir ../output/results_exhibits_xs/figures 
cap mkdir ../output/results_exhibits_xs/fig_data 


graph set window fontface default
graph set ps fontface default
graph set window fontfacemono "Consolas"
graph set ps fontfacemono "Consolas"
set scheme s1color


program main
	cf_scatter // JPE fig 8
	test_rates_figs // JPE fig 7
	vary_lb_cf // JPE fig 9
	comp_statics // JPE fig 6
	cap log close
end


program cf_scatter // JPE fig 8
	qui import delimited "$RESULTS/estimate_w_xs/cf_outcomes.csv", encoding(ISO-8859-2) clear
	
	rename v1 cf_name
	drop if cf_name == "Actual"
	gen num = _n 
	replace cf_name = "Out-of-pocket" if cf_name == "Out-of-pocket NIPT"
	* Insurance coverage of NIPT scatters
	gen ins_cf = 1 if inlist(cf_name, "Invasive only", "Free NIPT", "Out-of-pocket", ///
	  "[1:51-1:200]", " > 1:201")
	// twoway (scatter consumersurplus gcostpc if ins_cf == 1, msize(small)) ///
	// (scatter consumersurplus gcostpc if ins_cf == 1, msymbol(i) mlabel(cf_name) mlabpos(2)), ///
	// yscale(r(200 500)) ylab(#5, format(%5.0fc)) ///
	// xscale(r(100 600)) xlab(100(100)600) ///
	// ytitle("Avg. consumer surpulus") xtitle("Government cost per pregnancy") ///
	// legend(off) name(cf_scatter_ins_cs)
	twoway (scatter consumersurplus gcostpc if ins_cf == 1, msize(small)) ///
	(scatter consumersurplus gcostpc if ins_cf == 1, msymbol(i) mlabel(cf_name) mlabcolor(black) mlabpos(2)), ///
	ylab(#5, format(%5.0fc)) ///
	xscale(r(0 600)) xlab(0(100)600) ///
	ytitle("Avg. consumer surpulus") xtitle("Government cost per pregnancy") ///
	legend(off) name(cf_scatter_ins_cs)
 
	graph export "${RESULTS}/results_exhibits_xs/figures/cf_scatter_ins_cs.pdf", ///
	  name(cf_scatter_ins_cs) as(pdf) replace 
end

program test_rates_figs // JPE fig 7
	use "$TEMP/estimate_w_xs/inv_only_cf_out.dta", clear
	keep sim_wgt sim_nipt sim_invasive p_i a_i c_i oop_i
	rename sim_nipt nipt_inv_only
	rename sim_invasive invasive_inv_only
	foreach cf in "free_nipt" "oop" "g_200" {
		
		* THIS IS AN _N MERGE. DON'T SORT DATA DIFFERENTLY BETWEEN CF's in part5b-counter_rec.R
		merge 1:1 _n using "$TEMP/estimate_w_xs/`cf'_cf_out.dta", assert(match) nogen ///
		  keepusing(sim_nipt sim_invasive )
		rename sim_nipt nipt_`cf' 
		rename sim_invasive invasive_`cf'
		
	}
	
	qui gen fetus_risk = 1 / p_i
	qui gen bin_number = .
	local count = 1
	local bin_size = 25
	forval bin_max = `bin_size'(`bin_size')1000 {
		local bin_inf = `bin_max' - `bin_size'
		qui replace bin_number = `count' if fetus_risk <= `bin_max' & fetus_risk > `bin_inf'
		local count = `count' + 1
	}
	
	local ytext = 1
	
	collapse (mean) nipt_inv_only nipt_free_nipt nipt_oop nipt_g_200 ///
	  invasive_inv_only invasive_free_nipt invasive_oop invasive_g_200 [aweight = sim_wgt], ///
	  by(bin_number)
	qui replace bin_number = bin_number * 25
	
	twoway (connect nipt_free_nipt bin_number, m(t) mcolor(blue) lcolor(blue)) ///
	  (connect nipt_oop bin_number, m(s) mcolor(orange) lcolor(orange)), ///
	  legend(on) name(nipt_rate_3cfs) ///
	  xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
	  yscale(r(0 1)) ytitle("", size(small)) ///
	  ylabel(0(.1)1, format(%3.1f)) ///
	  xlabel(2 "1/2" 51 "1/51" 200 "1/200" 400 "1/400" 600 "1/600" ///
		800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
	  legend(order(1 "Free NIPT" 2 "NIPT Out-of-pocket") ///
	    pos(6) span row(2) col(2)) ///
	  xline(51 200) ///
	  plotregion(margin(zero)) xscale(reverse) ///
	  text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
	  text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
	  text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall))	
	graph export "${RESULTS}/results_exhibits_xs/figures/nipt_rate_3cfs.pdf", ///
	  name(nipt_rate_3cfs) as(pdf) replace

	twoway (connect invasive_free_nipt bin_number, m(t) mcolor(blue) lcolor(blue)) ///
	  (connect invasive_oop bin_number, m(s) mcolor(orange) lcolor(orange)), ///
	  legend(on) name(invasive_rate_3cfs) ///
	  xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
	  yscale(r(0 1)) ytitle("", size(small)) ///
	  ylabel(0(.10)1, format(%3.1f)) ///
	  xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" ///
		800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
	  legend(order(1 "Free NIPT" 2 "NIPT Out-of-pocket") ///
	    pos(6) span row(2) col(2)) ///
	  xline(51 200) ///
	  plotregion(margin(zero)) xscale(reverse) ///
	  text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
	  text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
	  text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall))	
	  
	graph export "${RESULTS}/results_exhibits_xs/figures/invasive_rate_3cfs.pdf", ///
	  name(invasive_rate_3cfs) as(pdf) replace
	
	gsort -bin_number
	export delimited using "$RESULTS/results_exhibits_xs/fig_data/test_rates_3cfs.csv", replace
end

program vary_lb_cf // JPE fig 9
	qui import delimited $RESULTS/estimate_w_xs/vary_lb_cf.csv, encoding(ISO-8859-2) clear
	drop v1 
	rename (v2 v3 v4 v5 v6 v7 v8 v9 v10 v11) (lb share_bad gcost subt nipt_only inv_only both_tests term_no live_ca_pref_a cons_surplus)
	
	foreach var in subt nipt_only inv_only both_tests term_no live_ca_pref_a share_bad {
		replace `var' = `var' / 100
	}
	
	* Get normalized cons_surplus and gcost
	qui sum cons_surplus if lb == 2
	local cs_norm = r(mean)
	qui sum gcost if lb == 2
	local g_norm = r(mean)
	gen cons_surplus_norm = cons_surplus - `cs_norm'
	gen gcost_norm = gcost - `g_norm'
	
	* LB that maximizes CS - G
	qui gen soc_welfare = cons_surplus - 1.3*gcost
	qui sort soc_welfare
	qui sum lb if _n == _N
	local lb_max_sw = r(mean)
	di "lb that maximizes CS - G is"
	di "`lb_max_sw'"
	
	// twoway scatter cons_surplus soc_welfare lb, yaxis(1) msize(small) m(o t)  ylabel(-400(100)400, format(%4.0f) axis(1)) ytitle("Surplus and Welfare", axis(1)) ///
	// 	xtitle("Lower bound risk for NIPT coverage", size(small)) xscale(r(0 1000) reverse) yscale(r(-400 400) axis(1)) ///
	// 	xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" 800 "1/800" 1000 "1/1000") ///
	// 	 || scatter gcost lb, yaxis(2) msize(small) m(d) ytitle("Cost", axis(2)) ylabel(0(100)600, format(%4.0f) axis(2)) ///
	// 	xtitle("Lower bound risk for NIPT coverage", size(small)) xscale(r(0 1000) reverse)  yscale(r(0 600) axis(2)) ///
	//   xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" 800 "1/800" 1000 "1/1000") ///
	//   xline(`lb_max_sw') ///
	//   ||, legend(order(1 "Avg. Consumer Surplus" 2 "Social Welfare (CS - 1.3G)" 3 "G cost p.c.") pos(6) span row(1) col(3)) /// 
	//   name(vary_lb_cf) scale(0.8)
	twoway scatter cons_surplus soc_welfare lb, yaxis(1) msize(small) m(o t)  ylabel(-400(400)2800, format(%4.0f) axis(1)) ytitle("Surplus and Welfare", axis(1)) ///
		xtitle("Lower bound risk for NIPT coverage", size(small)) xscale(r(0 1000) reverse) yscale(r(-400 1200) axis(1)) ///
		xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" 800 "1/800" 1000 "1/1000") ///
		 || scatter gcost lb, yaxis(2) msize(small) m(d) ytitle("Cost", axis(2)) ylabel(0(100)600, format(%4.0f) axis(2)) ///
		xtitle("Lower bound risk for NIPT coverage", size(small)) xscale(r(0 1000) reverse)  yscale(r(0 600) axis(2)) ///
	  xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" 800 "1/800" 1000 "1/1000") ///
	  xline(`lb_max_sw') ///
	  ||, legend(order(1 "Avg. Consumer Surplus" 2 "Social Welfare (CS - 1.3G)" 3 "G cost p.c.") pos(6) span row(1) col(3)) /// 
	  name(vary_lb_cf) scale(0.8)
	
	graph export "${RESULTS}/results_exhibits_xs/figures/vary_lb_cf.pdf", as(pdf) replace	
	
	twoway (lpoly live_ca_pref_a lb) (lpoly term_no lb), yscale(r(0 0.0022)) ///
	    ylabel(0(0.0005)0.002, format(%5.4fc)) ///
	 xtitle("Lower bound risk for NIPT coverage", size(small)) xscale(r(200 1000) reverse) ///
	  xlabel(50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" 800 "1/800" 1000 "1/1000") xline(200) ///
	  legend(order(1 "Live, chrom ab, a > c" 2 "Terminated, no chrom ab") ///
	    pos(6) span row(1) col(2)) ///
	  name(vary_lb_cf_out)
	graph export "${RESULTS}/results_exhibits_xs/figures/vary_lb_cf_out.pdf", as(pdf) replace
	 
	gsort -lb
	export delimited $RESULTS/results_exhibits_xs/vary_lb_cf_data.csv, replace 
end

program vary_ub_cf
	qui import delimited $RESULTS/estimate_w_xs/vary_ub_cf.csv, encoding(ISO-8859-2) clear
	drop v1 
	rename (v2 v3 v4 v5 v6 v7 v8 v9 v10) (ub share_bad gcost subt nipt_only inv_only both_tests term_no live_ca_pref_a )
	
	foreach var in subt nipt_only inv_only both_tests term_no live_ca_pref_a share_bad {
		replace `var' = `var' / 100
	}
	twoway (scatter share_bad ub, msize(small) m(d) yaxis(1) yscale(r(0 0.002) axis(1)) ///
	    ylabel(0(0.0005)0.002, format(%5.4fc) axis(1)) ///
	    ytitle("Share bad outcomes (%)", axis(1) orientation(vertical) size(small))) ///
	  (scatter gcost ub, msize(small) m(0) yaxis(2) yscale(r(0 300) axis(2)) ///
	    ylabel(0(50)300, format(%4.0f) axis(2)) ///
	    ytitle("Government cost per pregnancy", axis(2) orientation(rvertical) size(small))), ///
	  xtitle("Upper bound risk for NIPT coverage", size(small)) xscale(r(2 51) reverse) ///
	   xlabel(2 "1/2" 25 "1/25" 50 "1/50" 75 "1/75" 100 "1/100" 125 "1/125" 150 "1/150" 175 "1/175" ///
	     200 "1/200") xline(51) /// 
	  legend(order(1 "Share bad outcomes" 2 "G cost p.c.") pos(6) span row(1) col(2)) ///
	  name(vary_ub_cf)
	graph export "${RESULTS}/results_exhibits_xs/figures/vary_ub_cf.pdf", as(pdf) replace
	
	twoway scatter share_bad live_ca_pref_a term_no ub, msize(small) m(d o t) yscale(r(0 0.0025)) ///
	    ylabel(0(0.0005)0.0025, format(%5.4fc)) ///
	  ytitle("Share (%)",orientation(vertical) size(small)) ///
	 xtitle("Upper bound risk for NIPT coverage", size(small)) xscale(r(2 51) reverse) ///
	   xlabel(2 "1/2" 25 "1/25" 50 "1/50" 75 "1/75" 100 "1/100" 125 "1/125" 150 "1/150" 175 "1/175" ///
	     200 "1/200") xline(51) /// 
	  legend(order(1 "Bad outcome" 2 "Live, chrom ab, a > c" 3 "Terminated, no chrom ab") ///
	    pos(6) span row(1) col(3)) ///
	  name(vary_ub_cf_out)
	graph export "${RESULTS}/results_exhibits_xs/figures/vary_ub_cf_out.pdf", as(pdf) replace
	
	twoway scatter subt nipt_only inv_only both_tests ub, msize(small) m(o d t s) yscale(r(0 0.40)) ///
	  ylabel(0(.10).40, format(%3.2fc)) ///
	  ytitle("Share (%)",orientation(vertical) size(small)) ///
	 xtitle("Upper bound risk for NIPT coverage", size(small)) xscale(r(2 51) reverse) ///
	   xlabel(2 "1/2" 25 "1/25" 50 "1/50" 75 "1/75" 100 "1/100" 125 "1/125" 150 "1/150" 175 "1/175" ///
	     200 "1/200") xline(51) /// 
	  legend(order(1 "Any subsequent testing" 2 "NIPT only" 3 "Invasive only" 4 "NIPT followed by invasive") pos(6) span row(1) col(3)) ///
	  name(vary_ub_cf_test)
	graph export "${RESULTS}/results_exhibits_xs/figures/vary_ub_cf_test.pdf", as(pdf) replace
end

program comp_statics
	use "${TEMP}/results_exhibits_xs/comp_model.dta", clear
	qui gen fetus_risk = 1 / p_i
	qui gen bin_number = .
	local count = 1
	local bin_size = 25
	forval bin_max = `bin_size'(`bin_size')1000 {
		local bin_inf = `bin_max' - `bin_size'
		qui replace bin_number = `count' if fetus_risk <= `bin_max' & fetus_risk > `bin_inf'
		local count = `count' + 1
	}
	
	* From drawn a and c
	local ytext = 800
	preserve 
	collapse (mean) wtp_nipt, by(bin_number) 
	qui replace bin_number = bin_number * 25
	twoway connect wtp_nipt bin_number, ///
	  legend(on) name(avg_wtp_nipt) ///
	  xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
	  yscale(r(0 650)) ytitle("", size(small)) ///
	  ylabel(0(100)600, format(%3.0f)) ///
	  xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" ///
		800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
	  m(o) mcolor(blue) lcolor(blue) ///
	  legend(off) ///
	  xline(51 200) ///
	  plotregion(margin(zero)) xscale(reverse) ///
	  text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
	  text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
	  text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall))  
	
	graph export "${RESULTS}/results_exhibits_xs/figures/avg_wtp_nipt.pdf", as(pdf) replace
	restore
	
	// local ytext = 1
	// gen share_zero = 1 * (wtp_nipt > 0)
	// gen share_cost = 1 * (wtp_nipt > 567.5)
	
	// preserve 
	// collapse (mean) share_zero share_cost, by(bin_number)
	// qui replace bin_number = bin_number * 25
	// twoway (connect share_zero bin_number, m(o) mcolor(blue) lcolor(blue)) ///
	//   (connect share_cost bin_number, m(t) mcolor(green) lcolor(green)), ///
	//   legend(on) name(share_wtp_nipt) ///
	//   xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
	//   yscale(r(0 1)) ytitle("", size(small)) ///
	//   ylabel(#10, format(%3.2f)) ///
	//   xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" ///
	// 	800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
	//   legend(order(1 "WTP > 0" 2 "WTP > cost") ///
	// 		pos(6) span row(1) col(3)) ///
	//   xline(51 200) ///
	//   plotregion(margin(zero)) xscale(reverse) ///
	//   text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
	//   text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
	//   text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall)) 
	
	// graph export "${RESULTS}/results_exhibits_xs/figures/share_wtp_nipt.pdf", as(pdf) replace
	// restore
	
	* Comp stat in a and c
	foreach k in a c {
		local ytext = 1000
		preserve 
		collapse (mean) wtp_nipt_`k'_10 wtp_nipt_`k'_50 wtp_nipt_`k'_75 wtp_nipt_`k'_90, by(bin_number)
		qui replace bin_number = bin_number * 25
	
		  twoway (connect wtp_nipt_`k'_50 bin_number, m(o) mcolor(blue) lcolor(blue)) ///
		  (connect wtp_nipt_`k'_75 bin_number, m(t) mcolor(green) lcolor(green)) ///
		  (connect wtp_nipt_`k'_90 bin_number, m(s) mcolor(orange) lcolor(orange)), /// 
		  legend(on) name(avg_wtp_nipt_`k') ///
		  xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
		  yscale(r(0 650)) ytitle("", size(small)) ///
		  ylabel(0(100)600, format(%3.0f)) ///
		  xlabel(2 "1/2" 51 "1/51" 200 "1/200" 400 "1/400" 600 "1/600" ///
			800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
		  legend(order(1 "`k' - 50th percentile" 2 "`k' - 75th percentile" 3 "`k' - 90th percentile") ///
			pos(6) span row(1) col(3)) ///
		  xline(51 200) ///
		  plotregion(margin(zero)) xscale(reverse) ///
		  text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
		  text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
		  text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall))
		  
		graph export "${RESULTS}/results_exhibits_xs/figures/avg_wtp_nipt_`k'.pdf", as(pdf) replace
		restore
		
		// local ytext = 1
		// foreach lev in 10 50 75 90 {
		// 	gen share_zero_`k'_`lev' = 1 * (wtp_nipt_`k'_`lev' > 0)
		// 	gen share_cost_`k'_`lev'= 1 * (wtp_nipt_`k'_`lev' > 567.5)
		// }
		// preserve 
		// collapse (mean) share_zero_`k'_10 share_zero_`k'_50 share_zero_`k'_75 share_zero_`k'_90 ///
		//   share_cost_`k'_10 share_cost_`k'_50 share_cost_`k'_75 share_cost_`k'_90, by(bin_number)
		// qui replace bin_number = bin_number * 25
		// twoway scatter share_zero_`k'_10 share_zero_`k'_50 share_zero_`k'_75 share_zero_`k'_90 ///
		//   share_cost_`k'_10 share_cost_`k'_50 share_cost_`k'_75 share_cost_`k'_90 bin_number, ///
		//   legend(on) name(share_wtp_nipt_`k') ///
		//   xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
		//   yscale(r(0 1)) ytitle("", size(small)) ///
		//   ylabel(#10, format(%3.2f)) ///
		//   xlabel(2 "1/2" 51 "1/51" 200 "1/200" 400 "1/400" 600 "1/600" ///
		// 	800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
		//   m(o o o o t t t t) mcolor(blue green orange red blue green orange red) ///
		//   legend(order(1 "WTP > 0: `k' - 10th percentile" 2 "WTP > 0: `k' - 50th percentile" 3 "WTP > 0: `k' - 75th percentile" 4 "WTP > 0: `k' - 90th percentile" ///
		//     5 "WTP > cost: `k' - 10th percentile" 6 "WTP > cost: `k' - 50th percentile" 7 "WTP > cost: `k' - 75th percentile" 8 "WTP > cost: `k' - 90th percentile") ///
		// 	pos(6) span row(2) col(3)) ///
		//   xline(51 200) ///
		//   plotregion(margin(zero)) xscale(reverse) ///
		//   text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
		//   text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
		//   text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall))
		  
		// graph export "${RESULTS}/results_exhibits_xs/figures/share_wtp_nipt_`k'.pdf", as(pdf) replace
		// restore 
	}
	
	// * Comp stat in a and c splitting share zero and share cost
	// foreach k in a c {
	// 	di "`k'"
	// 	local ytext = 1
	// 	preserve
	// 	collapse (mean) share_zero_`k'_10 share_zero_`k'_50 share_zero_`k'_75 share_zero_`k'_90 ///
	// 	  share_cost_`k'_10 share_cost_`k'_50 share_cost_`k'_75 share_cost_`k'_90, by(bin_number)
	// 	qui replace bin_number = bin_number * 25
	// 	scatter share_zero_`k'_10 share_zero_`k'_50 share_zero_`k'_75 share_zero_`k'_90 bin_number, ///
	// 	  legend(on) name(share0_wtp_nipt_`k') ///
	// 	  xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
	// 	  yscale(r(0 1)) ytitle("", size(small)) ///
	// 	  ylabel(#10, format(%3.2f)) ///
	// 	  xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" ///
	// 		800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
	// 		 m(o t s d) mcolor(blue green orange red) ///
	// 	  legend(order(1 "`k' - 10th percentile" 2 "`k' - 50th percentile"   3 "`k' - 75th percentile" ///
	// 	    4 "`k' - 90th percentile") ///
	// 		pos(6) span row(1) col(3)) ///
	// 	  xline(51 200) ///
	// 	  plotregion(margin(zero)) xscale(reverse) ///
	// 	  text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
	// 	  text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
	// 	  text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall)) 
		  
	// 	graph export "${RESULTS}/results_exhibits_xs/figures/share0_wtp_nipt_`k'.pdf", ///
	// 	  name(share0_wtp_nipt_`k') as(pdf) replace
		
	// 	scatter share_cost_`k'_10 share_cost_`k'_50 share_cost_`k'_75 share_cost_`k'_90 bin_number, ///
	// 	  legend(on) name(shareCost_wtp_nipt_`k') ///
	// 	  xtitle("NT Screening result (estimated risk of chromosomal abnormalities)") ///
	// 	  yscale(r(0 1)) ytitle("", size(small)) ///
	// 	  ylabel(#10, format(%3.2f)) ///
	// 	  xlabel(2 "1/2" 50 "1/50" 200 "1/200" 400 "1/400" 600 "1/600" ///
	// 		800 "1/800" 1000 "1/1000", format(%7.0fc)) ///
	// 	  m(o t s d) mcolor(blue green orange red) ///
	// 	  legend(order(1 "`k' - 10th percentile" 2 "`k' - 50th percentile" 3 "`k' - 75th percentile" ///
	// 	    4 "`k' - 90th percentile") ///
	// 		pos(6) span row(1) col(3)) ///
	// 	  xline(51 200) ///
	// 	  plotregion(margin(zero)) xscale(reverse) ///
	// 	  text(`ytext' 25 "High risk", yaxis(1) place(n) size(vsmall)) ///
	// 	  text(`ytext' 125 "Medium risk", yaxis(1) place(n) size(vsmall)) ///
	// 	  text(`ytext' 600 "Low risk", yaxis(1) place(n) size(vsmall)) 
	// 	graph export "${RESULTS}/results_exhibits_xs/figures/shareCost_wtp_nipt_`k'.pdf", /// 
	// 		name(shareCost_wtp_nipt_`k') as(pdf) replace
	// 	restore 
	// }

	// collapse (mean) wtp_nipt* share_zero*, by(bin_number)
	// rename wtp_nipt avg_wtp_nipt
	// foreach k in a c {
	// 	foreach lev in 10 50 75 90 {
	// 		rename wtp_nipt_`k'_`lev' avg_wtp_nipt_`k'_`lev'
	// 	}
	// }
	// gsort -bin_number
	// qui export delimited $RESULTS/results_exhibits_xs/fig_data/wtp_nipt_fig_data.csv, replace 
	
	// *rm "${TEMP}/results_exhibits/comp_model.dta"
	
	

end


* EXECUTE
main
